
# Setup ------------------------------------------------------------------------
# Options
options(stringsAsFactors = FALSE)
# Packages
library(pacman)
# p_load(HiClimR, data.table, stringr, lubridate, lfe, fixest, magrittr, here)
p_load(data.table, stringr, lubridate, fixest, magrittr, here)
# Directory of joined bill-weather files
dir_project = "T:/Home/Edward/Elasticity"
dir_csv     = file.path(dir_project, "DataCsv")
dir_rds     = file.path(dir_project, "DataR")
dir_results = file.path(dir_rds, "Results/Results20200818/")

# Load the common zip code panel -----------------------------------------------
# Load the common zip code data
bill_dt = fst::read_fst(
  # file.path(dir_rds, "Prices/Joint/commonCitiesZipsBills.fst"),
  file.path(dir_rds, "Prices/Joint/zipsNeighbor0.fst"),
  as.data.table = TRUE
)
# # Load 'neighbor' datasets
# bill_dt = lapply(
#   X = c(
#     file.path(dir_rds, "Prices/Joint/zipsNeighbor0.fst"),
#     file.path(dir_rds, "Prices/Joint/zipsNeighbor1.fst"),
#     file.path(dir_rds, "Prices/Joint/zipsNeighbor2.fst"),
#     file.path(dir_rds, "Prices/Joint/zipsNeighbor3.fst")
#   ),
#   FUN = function(i) {
#     # Load the file
#     tmp = fst::read_fst(i, as.data.table = TRUE)
#     # Add indicator
#     tmp[, nbr := str_sub(i, -5, -5)]
#     # Return file
#     return(tmp)
#   }
# ) %>% rbindlist(use.names = TRUE, fill = TRUE)
# Drop unwanted variables
to_drop = setdiff(
  names(bill_dt),
  c(
    'begin_date', 'hh_id', 'nbr', 'thm_day', 'bill_hdd', 'city_ym',
    'price_mrg_lag2', 'spot_price_1week_lag2',
    'utility', 'care', 'season',
    'price_cluster',
    'flag_rev', 'flag_thm', 'bill_flag', 'days', 'days_lag1', 'days_lag2',
    'thm', 'thm_lag1', 'thm_lag2'
  )
)
bill_dt[, (to_drop) := NULL]

# Add leads and lags -----------------------------------------------------------
# Create lags and leads of flag_thm and flag_rev
setorder(bill_dt, hh_id, begin_date)
setkey(bill_dt, hh_id, begin_date)
bill_dt[, `:=`(
  # flag_rev_lead3 = shift(flag_rev,  3L, fill = NA, type = "lead"),
  # flag_thm_lead3 = shift(flag_thm,  3L, fill = NA, type = "lead"),
  # flag_bill_lead3 = shift(bill_flag,  3L, fill = NA, type = "lead"),
  # flag_rev_lead2 = shift(flag_rev,  2L, fill = NA, type = "lead"),
  # flag_thm_lead2 = shift(flag_thm,  2L, fill = NA, type = "lead"),
  # flag_bill_lead2 = shift(bill_flag,  2L, fill = NA, type = "lead"),
  flag_rev_lead1 = shift(flag_rev,  1L, fill = NA, type = "lead"),
  flag_thm_lead1 = shift(flag_thm,  1L, fill = NA, type = "lead"),
  flag_bill_lead1 = shift(bill_flag,  1L, fill = NA, type = "lead"),
  flag_rev_lag1  = shift(flag_rev,  1L, fill = NA, type = "lag"),
  flag_thm_lag1  = shift(flag_thm,  1L, fill = NA, type = "lag"),
  flag_bill_lag1  = shift(bill_flag,  1L, fill = NA, type = "lag"),
  flag_rev_lag2  = shift(flag_rev,  2L, fill = NA, type = "lag"),
  flag_thm_lag2  = shift(flag_thm,  2L, fill = NA, type = "lag"),
  flag_bill_lag2  = shift(bill_flag,  2L, fill = NA, type = "lag"),
  flag_rev_lag3  = shift(flag_rev,  3L, fill = NA, type = "lag"),
  flag_thm_lag3  = shift(flag_thm,  3L, fill = NA, type = "lag"),
  flag_bill_lag3  = shift(bill_flag,  3L, fill = NA, type = "lag")
  # flag_rev_lag4  = shift(flag_rev,  3L, fill = NA, type = "lag"),
  # flag_thm_lag4  = shift(flag_thm,  3L, fill = NA, type = "lag")
  # flag_bill_lag4  = shift(bill_flag,  3L, fill = NA, type = "lag")
), by = hh_id]
# Impose sample restrictions
sub_dt = bill_dt[(
  between(days, 28, 34) &
  between(days_lag1, 28, 34) &
  between(days_lag2, 28, 34) &
  thm > 0 & flag_thm == 0 & bill_flag == F & flag_rev == 0 &
  thm_lag1 > 0 & flag_thm_lag1 == 0 & flag_bill_lag1 == F & flag_rev_lag1 == 0 &
  thm_lag2 > 0 & flag_thm_lag2 == 0 & flag_bill_lag2 == F & flag_rev_lag2 == 0 &
  T
)]
rm(bill_dt); invisible(gc())

# # Impose sample restrictions ---------------------------------------------------
# # Drop observations flagged due to problems re-calculating bills
# # Keep bills with 28-34 days
# # NOTE: I flagged bills with absolute percent error of 5% or more.
# sub_dt = bill_dt[(
#   # between(days_lead1, 28, 34) &
#   between(days, 28, 34) &
#   between(days_lag1, 28, 34) &
#   between(days_lag2, 28, 34) &
#   # between(days_lag3, 28, 34) &
#   # thm_lead1 > 0 & flag_thm_lead1 == 0 & bill_flag_lead1 == F & flag_rev_lead1 == 0 &
#   thm > 0 & flag_thm == 0 & bill_flag == F & flag_rev == 0 &
#   thm_lag1 > 0 & flag_thm_lag1 == 0 & bill_flag_lag1 == F & flag_rev_lag1 == 0 &
#   thm_lag2 > 0 & flag_thm_lag2 == 0 & bill_flag_lag2 == F & flag_rev_lag2 == 0 &
#   # thm_lag3 > 0 & flag_thm_lag3 == 0 & bill_flag_lag3 == F & flag_rev_lag3 == 0 &
#   T
# )]
# rm(bill_dt); invisible(gc())

# Data work --------------------------------------------------------------------
# Add indicators for dimensions of heterogeneity
sub_dt[, `:=`(
  care = 1 * care,
  noncare = 1 * (care == FALSE),
  summer = 1 * (season == 'summer'),
  winter = 1 * (season == 'winter'),
  pge = 1 * (utility == 'pge'),
  socal = 1 * (utility == 'socal')
)]

# # Save dataset for Max ---------------------------------------------------------
#   # Save CSV with 'standard' variables
#   fwrite(
#     x = sub_dt[, .(
#       hh_id, begin_date, city_ym, price_cluster,
#       thm_day, bill_hdd, days,
#       price_mrg_lag2, price_avg_lag2, price_avg_mrg_lag2, price_base_lag2, mrg_sim_iv_10_14_lag2,
#       spot_price_1week_lag2, utility,
#       care, noncare, summer, winter
#     )],
#     file = file.path(dir_project, 'for-max.csv')
#   )

# Pooled estimates: City -------------------------------------------------------
# Marginal price
est_mrg_city = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Average price
est_avg_city = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Average marginal
est_avg_mrg_city = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Baseline price
est_base_city = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_base_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Simulated marginal price
est_sim_city = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(mrg_sim_iv_10_14_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save 
qs::qsave(
  x = list(est_mrg_city, est_avg_city, est_avg_mrg_city, est_base_city, est_sim_city),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-city.qs')
)
# Table of results
etable(
  est_mrg_city, est_avg_city, est_avg_mrg_city, est_base_city, est_sim_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-city.tex')
)
etable(
  est_mrg_city, est_avg_city, est_avg_mrg_city, est_base_city, est_sim_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  stage = 1,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-city-stage1.tex')
)
# Clean up
rm(est_mrg_city, est_avg_city, est_avg_mrg_city, est_base_city, est_sim_city)
gc()

# Pooled estimates: Zip --------------------------------------------------------
# Marginal price
est_mrg_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Average price
est_avg_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Average marginal
est_avg_mrg_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_avg_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Baseline price
est_base_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_base_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Simulated marginal price
est_sim_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(mrg_sim_iv_10_14_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save 
qs::qsave(
  x = list(est_mrg_zip, est_avg_zip, est_avg_mrg_zip, est_base_zip, est_sim_zip),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-zip.qs')
)
# Table of results
etable(
  est_mrg_zip, est_avg_zip, est_avg_mrg_zip, est_base_zip, est_sim_zip,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-zip.tex')
)
etable(
  est_mrg_zip, est_avg_zip, est_avg_mrg_zip, est_base_zip, est_sim_zip,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  stage = 1,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-zip-stage1.tex')
)
# Clean up
rm(est_mrg_zip, est_avg_zip, est_avg_mrg_zip, est_base_zip, est_sim_zip)
gc()

# CARE vs Non-CARE: City -------------------------------------------------------
# CARE
est_care_city = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Non-CARE
est_noncare_city = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care != TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
int_care_city = feols(
  log(thm_day) ~
    I(bill_hdd/days * noncare) + I(bill_hdd/days * care) | 
    hh_id ^ care + city_ym ^ care |
    I(log(price_mrg_lag2) * noncare) + I(log(price_mrg_lag2) * care) ~ 
    I(spot_price_1week_lag2 * noncare) +
    I(spot_price_1week_lag2 * noncare):utility +
    I(spot_price_1week_lag2 * care) +
    I(spot_price_1week_lag2 * care):utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_care_city, est_noncare_city, int_care_city),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'care-city.qs')
)
# Table of results
etable(
  est_care_city, est_noncare_city, int_care_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'care-city.tex')
)
etable(
  est_care_city, est_noncare_city, int_care_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  stage = 1,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'care-city-stage1.tex')
)
# Test coefficients for equality
# int_care_city$coefficients[1] - int_care_city$coefficients[2]
# sqrt(int_care_city$cov.scaled[1,1] + int_care_city$cov.scaled[2,2] - 2 * int_care_city$cov.scaled[1,2])
# Clean up
rm(est_care_city, est_noncare_city, int_care_city)
gc()

# CARE vs Non-CARE: City, Avg Price --------------------------------------------
# CARE
est_care_city = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Non-CARE
est_noncare_city = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care != TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
int_care_city = feols(
  log(thm_day) ~
    I(bill_hdd/days * noncare) + I(bill_hdd/days * care) | 
    hh_id ^ care + city_ym ^ care |
    I(log(price_avg_lag2) * noncare) + I(log(price_avg_lag2) * care) ~ 
    I(spot_price_1week_lag2 * noncare) +
    I(spot_price_1week_lag2 * noncare):utility +
    I(spot_price_1week_lag2 * care) +
    I(spot_price_1week_lag2 * care):utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_care_city, est_noncare_city, int_care_city),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'care-city-avg.qs')
)
# Table of results
etable(
  est_care_city, est_noncare_city, int_care_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'care-city-avg.tex')
)
etable(
  est_care_city, est_noncare_city, int_care_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  stage = 1,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'care-city-avg-stage1.tex')
)
# Test coefficients for equality
# int_care_city$coefficients[1] - int_care_city$coefficients[2]
# sqrt(int_care_city$cov.scaled[1,1] + int_care_city$cov.scaled[2,2] - 2 * int_care_city$cov.scaled[1,2])
# Clean up
rm(est_care_city, est_noncare_city, int_care_city)
gc()

# CARE vs Non-CARE: Zip --------------------------------------------------------
# CARE
est_care_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Non-CARE
est_noncare_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care != TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
int_care_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days * noncare) + I(bill_hdd/days * care) | 
    hh_id ^ care + zip_ym ^ care |
    I(log(price_mrg_lag2) * noncare) + I(log(price_mrg_lag2) * care) ~ 
    I(spot_price_1week_lag2 * noncare) +
    I(spot_price_1week_lag2 * noncare):utility +
    I(spot_price_1week_lag2 * care) +
    I(spot_price_1week_lag2 * care):utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_care_zip, est_noncare_zip, int_care_zip),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'care-zip.qs')
)
# Table of results
etable(
  est_care_zip, est_noncare_zip, int_care_zip,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'care-zip.tex')
)
etable(
  est_care_zip, est_noncare_zip, int_care_zip,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  stage = 1,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'care-zip-stage1.tex')
)
# Test coefficients for equality
# int_care_zip$coefficients[1] - int_care_zip$coefficients[2]
# sqrt(int_care_zip$cov.scaled[1,1] + int_care_zip$cov.scaled[2,2] - 2 * int_care_zip$cov.scaled[1,2])
# Clean up
rm(est_care_zip, est_noncare_zip, int_care_zip)
gc()

# Season: City -----------------------------------------------------------------
# Summer
est_summer_city =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[summer == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Winter
est_winter_city =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[summer != TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
int_season_city = feols(
  log(thm_day) ~
    I(bill_hdd/days * summer) + I(bill_hdd/days * winter) | 
    hh_id ^ season + city_ym ^ season |
    I(log(price_mrg_lag2) * summer) + I(log(price_mrg_lag2) * winter) ~ 
    I(spot_price_1week_lag2 * summer) +
    I(spot_price_1week_lag2 * summer):utility +
    I(spot_price_1week_lag2 * winter) +
    I(spot_price_1week_lag2 * winter):utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_summer_city, est_winter_city, int_season_city),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'season-city.qs')
)
# Table of results
etable(
  est_summer_city, est_winter_city, int_season_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'season-city.tex')
)
etable(
  est_summer_city, est_winter_city, int_season_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'season-city-stage1.tex')
)
# Test coefficients for equality
# int_season_city$coefficients[1] - int_season_city$coefficients[2]
# sqrt(int_season_city$cov.scaled[1,1] + int_season_city$cov.scaled[2,2] - 2 * int_season_city$cov.scaled[1,2])
# Clean up
rm(est_summer_city, est_winter_city, int_season_city)
gc()

# Season: City, average price --------------------------------------------------
# Summer
est_summer_city =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[summer == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Winter
est_winter_city =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[summer != TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
int_season_city = feols(
  log(thm_day) ~
    I(bill_hdd/days * summer) + I(bill_hdd/days * winter) | 
    hh_id ^ season + city_ym ^ season |
    I(log(price_avg_lag2) * summer) + I(log(price_avg_lag2) * winter) ~ 
    I(spot_price_1week_lag2 * summer) +
    I(spot_price_1week_lag2 * summer):utility +
    I(spot_price_1week_lag2 * winter) +
    I(spot_price_1week_lag2 * winter):utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_summer_city, est_winter_city, int_season_city),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'season-city-avg.qs')
)
# Table of results
etable(
  est_summer_city, est_winter_city, int_season_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'season-city-avg.tex')
)
etable(
  est_summer_city, est_winter_city, int_season_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'season-city-avg-stage1.tex')
)
# Test coefficients for equality
# int_season_city$coefficients[1] - int_season_city$coefficients[2]
# sqrt(int_season_city$cov.scaled[1,1] + int_season_city$cov.scaled[2,2] - 2 * int_season_city$cov.scaled[1,2])
# Clean up
rm(est_summer_city, est_winter_city, int_season_city)
gc()

# Season: Zip ------------------------------------------------------------------
# Summer
est_summer_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[summer == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Winter
est_winter_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[summer != TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
int_season_zip = feols(
  log(thm_day) ~
    I(bill_hdd/days * summer) + I(bill_hdd/days * winter) | 
    hh_id ^ season + zip_ym ^ season |
    I(log(price_mrg_lag2) * summer) + I(log(price_mrg_lag2) * winter) ~ 
    I(spot_price_1week_lag2 * summer) +
    I(spot_price_1week_lag2 * summer):utility +
    I(spot_price_1week_lag2 * winter) +
    I(spot_price_1week_lag2 * winter):utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_summer_zip, est_winter_zip, int_season_zip),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'season-zip.qs')
)
# Table of results
etable(
  est_summer_zip, est_winter_zip, int_season_zip,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'season-zip.tex')
)
etable(
  est_summer_zip, est_winter_zip, int_season_zip,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'season-zip-stage1.tex')
)
# Test coefficients for equality
# int_season_zip$coefficients[1] - int_season_zip$coefficients[2]
# sqrt(int_season_zip$cov.scaled[1,1] + int_season_zip$cov.scaled[2,2] - 2 * int_season_zip$cov.scaled[1,2])
# Clean up
rm(est_summer_zip, est_winter_zip, int_season_zip)
gc()

# Season-CARE interactions: City -----------------------------------------------
# Summer CARE
est_summer_care =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care == TRUE & summer == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Summer Non-CARE
est_summer_noncare =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[noncare == TRUE & summer == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Winter CARE
est_winter_care =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care == TRUE & winter == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Winter Non-CARE
est_winter_noncare =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[noncare == TRUE & winter == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
int_both = feols(
  log(thm_day) ~
    I(bill_hdd/days * summer * care) +
    I(bill_hdd/days * winter * care) +
    I(bill_hdd/days * summer * noncare) +
    I(bill_hdd/days * winter * noncare) | 
    hh_id ^ season ^ care + city_ym ^ season ^ care |
    I(log(price_mrg_lag2) * summer * care) +
    I(log(price_mrg_lag2) * winter * care) +
    I(log(price_mrg_lag2) * summer * noncare) +
    I(log(price_mrg_lag2) * winter * noncare) 
  ~ 
    I(spot_price_1week_lag2 * summer * care) +
    I(spot_price_1week_lag2 * summer * care):utility +
    I(spot_price_1week_lag2 * winter * care) +
    I(spot_price_1week_lag2 * winter * care):utility +
    I(spot_price_1week_lag2 * summer * noncare) +
    I(spot_price_1week_lag2 * summer * noncare):utility +
    I(spot_price_1week_lag2 * winter * noncare) +
    I(spot_price_1week_lag2 * winter * noncare):utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_summer_care, est_summer_noncare, est_winter_care, est_winter_noncare, int_both),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'both-city.qs')
)
# Table of results
etable(
  est_summer_care, est_summer_noncare, est_winter_care, est_winter_noncare, int_both,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'both-city.tex')
)
# Test Winter CARE vs Winter Non-CARE
# int_both$coefficients[2] - int_both$coefficients[4]
# sqrt(int_both$cov.scaled[2,2] + int_both$cov.scaled[4,4] - 2 * int_both$cov.scaled[2,4])
# Clean up
rm(est_summer_care, est_summer_noncare, est_winter_care, est_winter_noncare, int_both)
gc()

# Season-CARE interactions: City, Avg Price ------------------------------------
# Summer CARE
est_summer_care =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care == TRUE & summer == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Summer Non-CARE
est_summer_noncare =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[noncare == TRUE & summer == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Winter CARE
est_winter_care =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care == TRUE & winter == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Winter Non-CARE
est_winter_noncare =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[noncare == TRUE & winter == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
int_both = feols(
  log(thm_day) ~
    I(bill_hdd/days * summer * care) +
    I(bill_hdd/days * winter * care) +
    I(bill_hdd/days * summer * noncare) +
    I(bill_hdd/days * winter * noncare) | 
    hh_id ^ season ^ care + city_ym ^ season ^ care |
    I(log(price_avg_lag2) * summer * care) +
    I(log(price_avg_lag2) * winter * care) +
    I(log(price_avg_lag2) * summer * noncare) +
    I(log(price_avg_lag2) * winter * noncare) 
  ~ 
    I(spot_price_1week_lag2 * summer * care) +
    I(spot_price_1week_lag2 * summer * care):utility +
    I(spot_price_1week_lag2 * winter * care) +
    I(spot_price_1week_lag2 * winter * care):utility +
    I(spot_price_1week_lag2 * summer * noncare) +
    I(spot_price_1week_lag2 * summer * noncare):utility +
    I(spot_price_1week_lag2 * winter * noncare) +
    I(spot_price_1week_lag2 * winter * noncare):utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_summer_care, est_summer_noncare, est_winter_care, est_winter_noncare, int_both),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'both-city-avg.qs')
)
# Table of results
etable(
  est_summer_care, est_summer_noncare, est_winter_care, est_winter_noncare, int_both,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'both-city-avg.tex')
)
# Test Winter CARE vs Winter Non-CARE
# int_both$coefficients[2] - int_both$coefficients[4]
# sqrt(int_both$cov.scaled[2,2] + int_both$cov.scaled[4,4] - 2 * int_both$cov.scaled[2,4])
# Clean up
rm(est_summer_care, est_summer_noncare, est_winter_care, est_winter_noncare, int_both)
gc()

# Season-CARE interactions: Zip ------------------------------------------------
# Summer CARE
est_summer_care =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care == TRUE & summer == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Summer Non-CARE
est_summer_noncare =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[noncare == TRUE & summer == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Winter CARE
est_winter_care =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[care == TRUE & winter == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Winter Non-CARE
est_winter_noncare =  feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[noncare == TRUE & winter == TRUE],
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
int_both = feols(
  log(thm_day) ~
    I(bill_hdd/days * summer * care) +
    I(bill_hdd/days * winter * care) +
    I(bill_hdd/days * summer * noncare) +
    I(bill_hdd/days * winter * noncare) | 
    hh_id ^ season ^ care + zip_ym ^ season ^ care |
    I(log(price_mrg_lag2) * summer * care) +
    I(log(price_mrg_lag2) * winter * care) +
    I(log(price_mrg_lag2) * summer * noncare) +
    I(log(price_mrg_lag2) * winter * noncare) 
  ~ 
    I(spot_price_1week_lag2 * summer * care) +
    I(spot_price_1week_lag2 * summer * care):utility +
    I(spot_price_1week_lag2 * winter * care) +
    I(spot_price_1week_lag2 * winter * care):utility +
    I(spot_price_1week_lag2 * summer * noncare) +
    I(spot_price_1week_lag2 * summer * noncare):utility +
    I(spot_price_1week_lag2 * winter * noncare) +
    I(spot_price_1week_lag2 * winter * noncare):utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_summer_care, est_summer_noncare, est_winter_care, est_winter_noncare, int_both),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'both-zip.qs')
)
# Table of results
etable(
  est_summer_care, est_summer_noncare, est_winter_care, est_winter_noncare, int_both,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'both-zip.tex')
)
# Test Winter CARE vs Winter Non-CARE
# int_both$coefficients[2] - int_both$coefficients[4]
# sqrt(int_both$cov.scaled[2,2] + int_both$cov.scaled[4,4] - 2 * int_both$cov.scaled[2,4])
# Clean up
rm(est_summer_care, est_summer_noncare, est_winter_care, est_winter_noncare, int_both)
gc()

# Lags: City -------------------------------------------------------------------
# Lead 1
est_lead1 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lead1) ~ spot_price_1week_lead1 + spot_price_1week_lead1:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Lag 0 (contemporaneous)
est_lag0 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg) ~ spot_price_1week + spot_price_1week:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Lag 1
est_lag1 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag1) ~ spot_price_1week_lag1 + spot_price_1week_lag1:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Lag 2
est_lag2 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
all_lags = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lead1) +
    log(price_mrg) +
    log(price_mrg_lag1) +
    log(price_mrg_lag2) 
  ~ 
    spot_price_1week_lead1 + spot_price_1week_lead1:utility +
    spot_price_1week + spot_price_1week:utility +
    spot_price_1week_lag1 + spot_price_1week_lag1:utility +
    spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_lead1, est_lag0, est_lag1, est_lag2, all_lags),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'lags-city.qs')
)
# Table
etable(
  list(est_lead1, est_lag0, est_lag1, est_lag2, all_lags),
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'lags-city.tex')
)
# Clean up
rm(est_lead1, est_lag0, est_lag1, est_lag2, all_lags)
gc()

# Lags: City, Avg Price --------------------------------------------------------
# Lead 1
est_lead1 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lead1) ~ spot_price_1week_lead1 + spot_price_1week_lead1:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Lag 0 (contemporaneous)
est_lag0 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg) ~ spot_price_1week + spot_price_1week:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Lag 1
est_lag1 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag1) ~ spot_price_1week_lag1 + spot_price_1week_lag1:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Lag 2
est_lag2 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
all_lags = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_avg_lead1) +
    log(price_avg) +
    log(price_avg_lag1) +
    log(price_avg_lag2) 
  ~ 
    spot_price_1week_lead1 + spot_price_1week_lead1:utility +
    spot_price_1week + spot_price_1week:utility +
    spot_price_1week_lag1 + spot_price_1week_lag1:utility +
    spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_lead1, est_lag0, est_lag1, est_lag2, all_lags),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'lags-city-avg.qs')
)
# Table
etable(
  list(est_lead1, est_lag0, est_lag1, est_lag2, all_lags),
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'lags-city-avg.tex')
)
# Clean up
rm(est_lead1, est_lag0, est_lag1, est_lag2, all_lags)
gc()

# Lags: Zip --------------------------------------------------------------------
# Lead 1
est_lead1 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lead1) ~ spot_price_1week_lead1 + spot_price_1week_lead1:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Lag 0 (contemporaneous)
est_lag0 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg) ~ spot_price_1week + spot_price_1week:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Lag 1
est_lag1 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag1) ~ spot_price_1week_lag1 + spot_price_1week_lag1:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Lag 2
est_lag2 = est_mrg
# Lag 3
est_lag3 = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lag3) ~ spot_price_1week_lag3 + spot_price_1week_lag3:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt[(
    between(days_lag3, 28, 34) &
    thm_lag3 > 0 & flag_thm_lag3 == 0 & bill_flag_lag3 == F & flag_rev_lag3 == 0
  )],
  mem.clean = TRUE,
  lean = FALSE
)
# Joint estimation
all_lags = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + zip_ym |
    log(price_mrg_lead1) +
    log(price_mrg) +
    log(price_mrg_lag1) +
    log(price_mrg_lag2) +
    log(price_mrg_lag3)
  ~ 
    spot_price_1week_lead1 + spot_price_1week_lead1:utility +
    spot_price_1week + spot_price_1week:utility +
    spot_price_1week_lag1 + spot_price_1week_lag1:utility +
    spot_price_1week_lag2 + spot_price_1week_lag2:utility +
    spot_price_1week_lag3 + spot_price_1week_lag3:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save
qs::qsave(
  x = list(est_lead1, est_lag0, est_lag1, est_lag2, all_lags),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'lags-zip.qs')
)
# Table
etable(
  list(est_lead1, est_lag0, est_lag1, est_lag2, est_lag3, all_lags),
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'lags-zip.tex')
)
# Clean up
rm(est_lead1, est_lag0, est_lag1, est_lag2, est_lag3, all_lags)
gc()


# Pooled estimates: City, No HDDs ----------------------------------------------
# Marginal price
est_mrg_city = feols(
  log(thm_day) ~
    1 | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Average price
est_avg_city = feols(
  log(thm_day) ~
    1 | 
    hh_id + city_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Average marginal
est_avg_mrg_city = feols(
  log(thm_day) ~
    1 | 
    hh_id + city_ym |
    log(price_avg_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Baseline price
est_base_city = feols(
  log(thm_day) ~
    1 | 
    hh_id + city_ym |
    log(price_base_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Simulated marginal price
est_sim_city = feols(
  log(thm_day) ~
    1 | 
    hh_id + city_ym |
    log(mrg_sim_iv_10_14_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save 
qs::qsave(
  x = list(est_mrg_city, est_avg_city, est_avg_mrg_city, est_base_city, est_sim_city),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-city-nohdd.qs')
)
# Table of results
etable(
  est_mrg_city, est_avg_city, est_avg_mrg_city, est_base_city, est_sim_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-city-nohdd.tex')
)
etable(
  est_mrg_city, est_avg_city, est_avg_mrg_city, est_base_city, est_sim_city,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  stage = 1,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-city-nohdd-stage1.tex')
)
# Clean up
rm(est_mrg_city, est_avg_city, est_avg_mrg_city, est_base_city, est_sim_city)
gc()

# Pooled estimates: Zip, No HDDs -----------------------------------------------
# Marginal price
est_mrg_zip = feols(
  log(thm_day) ~
    1 | 
    hh_id + zip_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Average price
est_avg_zip = feols(
  log(thm_day) ~
    1 | 
    hh_id + zip_ym |
    log(price_avg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Average marginal
est_avg_mrg_zip = feols(
  log(thm_day) ~
    1 | 
    hh_id + zip_ym |
    log(price_avg_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Baseline price
est_base_zip = feols(
  log(thm_day) ~
    1 | 
    hh_id + zip_ym |
    log(price_base_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Simulated marginal price
est_sim_zip = feols(
  log(thm_day) ~
    1 | 
    hh_id + zip_ym |
    log(mrg_sim_iv_10_14_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
# Save 
qs::qsave(
  x = list(est_mrg_zip, est_avg_zip, est_avg_mrg_zip, est_base_zip, est_sim_zip),
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-zip-nohdd.qs')
)
# Table of results
etable(
  est_mrg_zip, est_avg_zip, est_avg_mrg_zip, est_base_zip, est_sim_zip,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-zip-nohdd.tex')
)
etable(
  est_mrg_zip, est_avg_zip, est_avg_mrg_zip, est_base_zip, est_sim_zip,
  style.tex = style.tex('aer'),
  fitstat = ~ r2 + n, tex = TRUE,
  stage = 1,
  file = file.path(dir_project, 'DataR', 'Results', 'Results20220818', 'pooled-zip-nohdd-stage1.tex')
)
# Clean up
rm(est_mrg_zip, est_avg_zip, est_avg_mrg_zip, est_base_zip, est_sim_zip)
gc()
